# Survival-analysis simulation (with censoring)
rm(list = ls())
library(Mediana)

# Define data model ----

# Endpoint parameters
hazard_ratio <- 0.72
OS.control <- 10 # Months
OS.test <- OS.control / hazard_ratio

hrate.control <- log(2) / OS.control # Hazard rate
hrate.test <- log(2) / OS.test

outcome.control <- parameters(rate = hrate.control) # Endpoint object 
outcome.test <- parameters(rate = hrate.test)

# Number of events parameters
event.count <- c(200, 340) # Total event count between 200 and 340
rand.ratio <- c(1, 1) # 1:1, 1:2 is c(1, 2)

# Dropout
dropoutrate <- parameters(rate = 0.015)


# Data model
dmTrial <- DataModel() +
  OutcomeDist(outcome.dist = "ExpoDist") +
  Event(n.events = event.count, rando.ratio = rand.ratio) +
  Sample(id = "Control",
         outcome.par = parameters(outcome.control)) +
  Sample(id = "Test",
         outcome.par = parameters(outcome.test)) +
  Design(enroll.period = 12,
         study.duration = 24,
         enroll.dist = "UniformDist",
         dropout.dist = "ExpoDist",
         dropout.dist.par = dropoutrate)


# Define analysis model ----
amTrial <- AnalysisModel() +
  Test(id = "Ctrl vs Test",
       samples = samples("Control", "Test"),
       method = "LogrankTest")

amTrial <- amTrial +
  Statistic(id = "Hazard ratio",
            samples = samples("Control", "Test"),
            method = "HazardRatioStat",
            par = parameters(method = "Cox")) +
  Statistic(id = "Events - Control",
            samples = samples("Control"),
            method = "EventCountStat") +
  Statistic(id = "Events - Test",
            samples = samples("Test"),
            method = "EventCountStat") +
  Statistic(id = "Subjects - Control",
            samples = samples("Control"),
            method = "PatientCountStat") +
  Statistic(id = "Subjects - Test",
            samples = samples("Test"),
            method = "PatientCountStat")

# Define evaluation model ----
emTrial <- EvaluationModel() +
  Criterion(id = "PowerAnalysis",
            method = "MarginalPower",
            tests = tests("Ctrl vs Test"),
            labels = c("Control vs test"),
            par = parameters(alpha = 0.05)) +
  Criterion(id = "HazardRatio",
            method = "MeanSumm",
            statistics = statistics("Hazard ratio"),
            labels = c("Average hazard ratio")) +
  Criterion(id = "Events",
            method = "MeanSumm",
            statistics = statistics("Events - Control", "Events - Test"),
            labels = c("Events (Control)", "Events (Test)")) +
  Criterion(id = "Events",
            method = "MeanSumm",
            statistics = statistics("Subjects - Control", "Subjects - Test"),
            labels = c("Events (Control)", "Events (Test)"))
  
            


# CSE ----
simparam <- SimParameters(n.sims = 10, proc.load = "full", seed = 1)
cseTrial <- CSE(dmTrial, amTrial, emTrial, simparam)

summary(cseTrial)
